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We study chaotic orbits of conservative low-dimensional maps and present numerical results 
showing that the probability density functions (pdfs) of the sum of N iterates in the large N limit 
exhibit very interesting time-evolving statistics. In some cases where the chaotic layers are thin 
and the (positive) maximal Lyapunov exponent is small, long-lasting quasi-stationary states (QSS) 
are found, whose pdfs appear to converge to g-Gaussians associated with nonextensive statistical 
mechanics. More generally, however, as N increases, the pdfs describe a sequence of QSS that pass 
from a ^-Gaussian to an exponential shape and ultimately tend to a true Gaussian, as orbits diffuse 
to larger chaotic domains and the phase space dynamics becomes more uniformly ergodic. 



I. INTRODUCTION 

As is well-known, invariant closed curves of area-preserving maps present complete barriers to orbits evolving 
inside resonance islands in the two-dimensional phase space. Outside these regions, there exist families of smaller 
islands and invariant Cantor sets (often called cantori), to which chaotic orbits are observed to "stick" for very long 
times. Thus, at the boundaries of these islands, an 'edge of chaos' develops with vanishing (or very small) Lyapunov 
exponents, where trajectories yield quasi-stationary states (QSS) that are often very long-lived. Such phenomena 
have been thoroughly studied to date in terms of a number of dynamical mechanisms responsible for chaotic transport 
in area-preserving maps and low-dimensional Hamiltonian systems (see e.g. jl4 l3l|). 

In this paper we study numerically the probability density functions (pdfs) of sums of iterates of QSS characterized 
by non- vanishing Lyapunov exponents, aiming to understand the connection between their intricate phase space 
dynamics and their time-evolving statistics. Our approach, therefore, is in the context of the Central Limit Theorem 
(CLT), although in many cases our pdfs do not converge to a single shape but pass through several ones. One case 
where convergence is known to exist is when the dynamics is bounded and uniformly hyperbolic (as e.g. in the case 
of Sinai billiards) and the associated pdf is a Gaussian. However, even in nonhyperbolic conservative models, there 
are regions where trajectories are essentially ergodic and mixing, so that Gaussians are ultimately observed, as the 
number of iterations grows. In such cases the maximal Lyapunov exponent is positive and bounded away from zero. 
What happens, however, when the motion is "weakly" chaotic and explores domains with intricate invariant sets, 
where the maximal (positive) Lyapunov exponent is very small? It is the purpose of this work to explore the statistics 
of such regions and determine the type of QSS generated by their dynamics. 

Recently, there has been a number of interesting studies of such pdfs of one-dimensional maps p], HO, IH, l24j 
and higher-dimensional conservative maps [18| in precisely 'edge of chaos' domains, where the maximal Lyapunov 
exponent either vanishes or is very close to zero. These studies provide evidence for the existence of q- Gaussian 
distributions, in the context of the Central Limit Theorem. This sparked off fierce controversy [ll| but, for one- 
dimensional maps, the argument has been resolved. In fact, 0, H3 undoubtedly show that, when approaching the 
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critical point while taking into account a proper scaling relation that involves the vicinity of the critical point and the 
Feigenbaum constant 5, the pdfs of sums of iterates of the logistic map are approximated by a q- Gaussian far better 
that the Levy distribution suggested in [11]. This suggests the need for a more thorough investigation of these systems 
within a nonextensive statistical mechanics approach, based on the nonadditive entropy S q |25l l26j . According to 
this approach, the pdfs optimizing (under appropriate constraints) S q are g-Gaussian distributions that represent 
metastable states @, HI, EBt EsJ , or QSS of the dynamics. 

The validity of a Central Limit Theorem (CLT) has been verified for deterministic systems 0, [H, [HI and, more 
recently, a ^-generalization of the CLT was proved demonstrating that, for certain classes of strongl y c orrelated 
random variables, their rescaled sums approach not a Gaussian, but a g-Gaussian limit distribution jl2[ |28[ |29| . 
Systems statistically described by power-law probability distributions (a special case of which are g-Gaussians) are in 
fact so ubiquitous [2l|, that some authors claimed that the normalization technique of a type of data that characterizes 
the measurement device is one of the reasons of their occurrence [30|: This is the case of normalized and centered 
sums of data that exhibit elliptical symmetry, but not necessarily the case of the iterates of deterministic maps, as 
can be inferred by the verification of a classical CLT for the paradigmatic example of the fully chaotic logistic map. 

In this paper, we follow this reasoning and compute first, in weakly chaotic domains of conservative maps, the pdf 
of the rescaled sum of N iterates, in the large N limit, and for many different initial conditions. We then connect our 
results with specific properties of the phase space dynamics of the maps and distinguish cases where the pdfs represent 
long-lived QSS described by g-Gaussians. We generally find that, as N grows, these pdfs pass from a g-Gaussian 
to an exponential (having a triangular shape in our semi- log plots), ultimately tending to become true Gaussians, as 
"stickiness" to cantori apparently subsides in favor of more uniformly chaotic (or ergodic) motion. 

In section II we begin our study by a detailed study of QSS, their pdfs and corresponding dynamics in two- 
dimensional Ikeda and MacMillan maps. In section III we briefly discuss analogous phenomena in 4-dimensional 
conservative maps and end with our conclusions in section IV. 



II. TWO-DIMENSIONAL AREA-PRESERVING MAPS 

Let us consider two-dimensional maps of the form: 

Vn+i = 9(x n ,y n ) 

and treat their chaotic orbits as generators of random variables. Even though this is not true for the iterates of a 
single orbit, we may still regard as random sequences those produced by many independently chosen initial conditions. 
In [l5|, the well known CLT assumption about the independence of N identically distributed random variables was 
replaced by a weaker property that essentially means asymptotic statistical independence. Thus, we may proceed to 
compute the generalized rescaled sums of their iterates Xi in the context of the classical CLT (see ((j ITHj ) : 

N 

Z N = N~^J2( x i-( x )) ( 2 ) 
i=i 

where (• • • ) implies averaging over a large number of iterations N and a large number of randomly chosen initial 
conditions Ni c . Due to the possible nonergodic and nonmixing behavior, averaging over initial conditions is an 
important ingredient of our approach. 

For fully chaotic systems (7 = 1/2), the distribution of ([2]) in the limit (TV — >> 00) is expected to be a Gaussian [l5| . 
Alternatively, however, we may define the non-rescaled variable zn 

N 

Z N = Y1 & ~ ( 3 ) 

i=l 

and analyze the probability density function (pdf) of z n normalized by its variance (so as to absorb the rescaling 
factor TV 7 ) as follows: 

First, we construct the sums Sjp obtained from the addition of TV x- iterates Xi (i = 0, . . . , TV) of the map (pQ) 

N 



stf = j24 j) (4) 



i=0 
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TABLE I: Maximal Lyapunov exponents of the Ikeda map, with C\ — 0.4, C2 = 6, R = 1 and u = 0.7, 0.8, 0.9, 1.0. 



u 


0.7 


0.8 


0.9 


1.0 


Lraax 


0.334 


0.344 


0.5076 


0.118 



where (j) represents the dependence of Sjp on the randomly chosen initial conditions Xq \ with j = 1, 2, N ic . Next, 
we focus on the centered and rescaled sums 



,0) 

N 




(s$ - (Sjp)) /*s = £ *« W "jEE ^ A™ ( 5 ) 



where a at is the standard deviation of the 5^ 



%C 3=1 

Next, we estimate the pdf of plotting the histograms of P(s$) f° r sufficiently small increments As^\= 0.05 
is used in all cases), so as to smoothen out fine details and check if they are well fitted by a g-Gaussian: 



P 



(s^)=P(0)[l^P(q-l)(s^)Y~ 9 (7) 



where q is the index of the nonadditive entropy S q and /3 is a 'inverse temperature' parameter. Note that as q — > 1 this 

distribution tends to a Gaussian, i.e., lim^i P(s$) = P(0)e _/3 ( s ^ ) 2 . From now on, we write z/a = s$ . We also 
remark that, due to the projection of the higher dimensional motion onto a single axis, the support of our distributions 
appears to consist of a dense set of values in z/a, although we cannot analytically establish its continuum nature. 



A. The Ikeda map 

Let us first examine by this approach the well-known Ikeda map Q: 

x n+1 = R + u(x n cos r - y n sin r) 
Vn+i = u(x n sin r + y n cos r) 

where r = C\ — Cij{\ + % 2 n + 2/n), R, u, C\, C2 are free parameters, and the Jacobian is J(R,u,r) = u 2 , so that (|8j) 
is dissipative for u < 1 and area-preserving for u = 1. This map was proposed as a model, under some simplifying 
assumptions, of the type of cell that might be used in an optical computer [2]. Fixing the values of C\ — 0.4, C2 = 6 
and R = 1 we observe that when u = 0.7,0.8,0.9, areas of the phase plane contract and strange attractors appear. 
In Fig. [T]we plot two different structures of the phase space dynamics for representative values of the parameter, u. 
The values of the positive (largest) Lyapunov exponent L max in these cases are listed in the Table IB 

Fig. [2] shows the corresponding pdf of the normalized variables (0 obtained for the two values of the parameter, 
u = 0.9, 1, in the large TV limit. In fact, we observe that for u = 0.7,0.8,0.9, the system possesses strange chaotic 
attractors whose pdfs can be properly fitted by Gaussians. These numerical results are not in disagreement with 
those of [22j| , on the 2-dimensional Henon map, where it was shown that its strange attr actor exhibits nonextensive 
properties (i.e., q ^ 1). In a fully chaotic domain, non-extensive properties need not be present and consequently 
pdfs of the sum of iterates should be Gaussian distributions. Now, for u = 0.7,0.8,0.9, the Ikeda map (j5J) generates 
strange attractors whose maximum Lyapunov exponent is positive and bounded away from zero (see Table [T|). This 
means that the motion is not at the 'edge of chaos' but rather in a chaotic sea and consequently the concepts involved 
in Boltzmann-Gibbs statistics are expected to hold. On the contrary, in the area-preserving case u = 1, the pdf of 
the sums of ([5|) converges to a non-Gaussian function (see Fig. Eb)- 

Now, in an 'edge of chaos' regime, one might expect to obtain a g-Gaussian limit distribution (q < 3), which 
generalizes Gaussians and extremizes the nonadditive entropy S q under appropriate constraints. Of course, the 
chaotic annulus shown in Fig. [1] for u = 1 does not represent an 'edge of chaos' regime, as the maximal Lyapunov 
exponent does not vanish (see Table [Ij) and the orbit appears to explore this annulus more or less uniformly. Hence 
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FIG. 1: Phase space plots of the Ikeda map for C\ — 0.4, C2 = 6, R = 1, and representative values of u. When u = 0.9, areas 
of the phase plane contract and a strange attractors appears. When u = 1, the map is area-preserving and a chaotic annular 
region is observed surrounding a domain about the origin where the motion is predominantly quasiperiodic. We use randomly 
chosen initial conditions from a square [0, 10 -4 ] x [0, 10 -4 ] about the origin (0,0). 




FIG. 2: (Color online) Pdfs of the normalized sums of iterates of the Ikeda map, for Ci = 0.4, C2 = 6, R = 1. N represents the 
number of (summed) iterates. Panel (a): Ni C is the number of randomly chosen initial conditions from the basin of attraction 
(dissipative case); black line corresponds to Gaussian function e - ^ 2 /^ , /3 — 0.5. Panel(b): Ni C is the number of randomly 
chosen initial conditions from a square [0, 1] x [0, 1] located inside the chaotic annular region of the area-preserving map; black 
line corresponds to (q = 5.3)-Gaussian functional. 



a g-Gaussian distribution in that case would not be expected. But appearances can be deceiving. The result we 
obtain is remarkable, as the central part of our pdf is well-fitted by a q- Gaussian functional with q = 5.3 up to very 
large N (see Fig JJJd). Although this is not a normalizable ^-Gaussian function (since q > 3 [26]), it is nevertheless 
striking enough to suggest that: (a) the motion within the annular region is not as uniformly ergodic as one might 
have expected and (b) the L max is not large enough to completely preclude 'edge of chaos' dynamics. 

All this motivated us to investigate more carefully similar phenomena in another family of area-preserving maps. 



B. The MacMillan map 



Consider the so-called perturbed MacMillan map, which may be interpreted as describing the effect of a simple 
linear focusing system supplemented by a periodic sequence of thin nonlinear lenses [17j: 

+ 2/i T ^ + e(^ n + /3x n ) (9j 
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FIG. 3: (Color online) Dynamical and statistical behavior of chaotic orbits of the MacMillan map for parameter values fi — 1.6, 
and e = 0.2, 0.9, 1.8 (from left to right). Figs.(a)-(c) represent the pdfs of the normalized sums of N iterates; Ni C is the number 
of randomly chosen initial conditions, from a square (0, 10 -6 ) x (0, 10 -6 ). Figs.(d)-(f) depict the corresponding phase space 
plots. 



TABLE II: Maximal Lyapunov exponents of the MacMillan map, with fi = 0.6 and e = 0.2, 0.5, 0.9, 1.1, 1.2, 1.8. 



e 


0.2 


0.5 


0.9 


1.1 


1.2 


1.8 


Lmax 


0.0867 


0.082 


0.0875 


0.03446 


0.0513 


0.05876 



where e, /?, /i are physically important parameters. The Jacobian is J(e, f3) = 1 — e/3, so that ((9]) is area-preserving for 
e/3 = 0, and dissipative for e/3 > 0. Here, we only consider the area-preserving case j3 — 0, so that the only relevant 
parameters are (e, fi). 

The unperturbed map yields a lemniscate invariant curve with a self-intersection at the origin that is a fixed point 
of saddle type. For e ^ 0, separatrices split and the map presents a thin chaotic layer around two islands. Increasing 
e, chaotic regions spread in the x ni y n plane. 

Within these chaotic regions, we have analyzed the histogram of the normalized sums of ([5]) for a wide range of 
parameters (e, fi) and have identified some generic pdfs in the form of q-Gaussians, and exponentials ~ e~ k ^ having 
a triangular shape on semi- logarithmic scale, which we call for convenience triangular distributions. Monitoring their 
'time evolution' under increasingly large numbers of iterations A/", we typically observe the occurrence of different 
QSS described by these distributions. We have also computed their L max and corresponding phase space plots and 
summarized our main results in Figures [3] and |4j The maximal Lyapunov exponents for the cases shown in Figures [3] 
and |U are listed in Table HI 

Below, we discuss the time-evolving statistics of two examples of the Mac Millan map, which represent respectively: 
(1) One set of cases with a 'figure eight' chaotic domain whose distributions pass through a succession of pfds before 
converging to an ordinary Gaussian (Fig. [3]), and (2) a set with more complicated chaotic domains extending around 
many islands, where q- Gaussian pdfs dominate the statistics for very long times and convergence to a Gaussian is not 
observed (Fig. |4]). 
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FIG. 4: (Color online) Dynamical and statistical behavior of chaotic orbits of the MacMillan map for parameter values fi — 1.6, 
and e = 0.5, 1.1, 1.2 (from left to right), where the orbits evolve around a central 'figure eight' chaotic region. Figs, (a)-(c) 
represent the pdfs of the normalized sums of N iterates; Ni C is the number of randomly chosen initial conditions, from a square 
(0, 10 -6 ) x (0, 10 -6 ). Figs, (d)-(f) depict the corresponding phase space plots. 



1. ( e = 0.9, fi - 1.6) -MacMillan map 

The (0.9, 1.6)-MacMillan map is a typical example producing time-evolving pdfs. As shown in Fig. [3j the corre- 
sponding phase space plots yield a seemingly simple chaotic region in the form of a 'figure eight' around two islands, yet 
the corresponding pdfs do not converge to a single distribution; rather they pass from a g-Gaussian-looking function 
to a triangular distribution. 

Analyzing carefully this time evolution of pdfs, we observed that there exist at least three long-lived QSS, whose 
iterates mix in the 2-dimensional phase space to generate superimposed pdfs of the corresponding sums ([5]). Conse- 
quently, for i = 1 . . . N = 2 16 , a QSS is produced whose pdf is close to a pure (q = 1.6)-Gaussian whose j3 parameter 
increases as N increases and the density of phase space plot grows (see Fig. EJ. This kind of distribution, in a fully 
chaotic region, is affected not only by a Lyapunov exponent being close to zero, but also by a "stickiness" effect around 
islands of regular motion. In fact, the boundaries of these islands is where the 'edge of chaos' regime is expected to 
occur in conservative maps (32j . 

Fig. [5] and Fig. [6] show some phase space plots for different numbers of iterates TV. Note that for N = 1 . . . 2 16 , 
these plots depict first a 'figure eight' chaotic region that evolves essentially around two islands (Fig.[5j). However, for 
N > 2 16 , a more complex structure emerges: Iterates stick around new islands, and an alternation of QSS is evident 
from g-Gaussian to exponentially decaying shapes (see Fig. [6]). 

Clearly, therefore, for e = 0.9 (and other similar cases with e = 0.2, 1.8) more than one QSS coexist whose pdfs 
are the superposition of their corresponding (q ^ 1)-Gaussians. Note in Fig. that this superposition of QSS occurs 
for 10 18 < N < 2 21 and produces a mixed distribution where the central part is still well-described by a (q = 1.6)- 
Gaussian. However, as we continue to iterate the map to N = 2 23 , this g-Gaussian is hidden by a superposition of 
intermediate states, passing to a triangular distribution. From here on, as N > 2 23 , the central part of the pdfs is 
close to a Gaussian (see Fig. [8] and Fig. [9j) and a true Gaussian is expected in the limit (TV — > oo). The evolution of 
this sequence of successive QSS as TV increases is shown in detail in Fig. [9] 
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FIG. 5: (Color online) Panel (a): PDFs of the renormalized sums of N iterates of the (e = 0.9, \i — 1.6)-MacMillan map, for 
N < 10 16 , and N% c randomly chosen initial condition in a square (0, 10 -6 ) x (0, 10 -6 ). Panel (bl)-(b2): Corresponding phase 
space plots for N — 2 12 and N — 2 16 . 
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FIG. 6: (e = 0.9, fi = 1.6)-MacMillan map partial phase space evolution. The iterates are calculated starting form a randomly 
chosen initial condition in a square (0, 10 -6 ) x (0, 10 -6 ). N is the number of plotted iterates. Note the long-standing quasi- 
stationary states that sequentially superimpose on phase space plots. 
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FIG. 7: Panel (a): (e = 0.9, \i — 1.6)-MacMillan phase space plots for i = 1 . . . N (N > 2 23 ) iterates, starting form a randomly 
chosen initial condition in a square (0, 10 -6 ) x (0, 10 -6 ). Panel (b)-(c): (Color online) Corresponding PDFs. Ni C is the number 
of randomly chosen initial condition in a square (0, 10 -6 ) x (0, 10 -6 ). 



2. (e = 1.2, /i = 1.6)-MacMillan map 



Let us now carefully analyze the behavior of the (1.2, 1.6)-MacMillan map, whose maximum Lyapunov exponent 
is L max « 0.05, smaller than that of the e = 0.9 case (L max w 0.08). As is clearly seen in Fig. [lOj a diffusive behavior 
sets in here that extends outward in phase space, envelopping a chain of islands of an order 8 resonance to which the 
orbits "stick" as the number of iterations grows to N = 2 19 . 

Again, chaotic motion starts by enveloping the same 'figure eight' as in the e = 0.9 case and the central part of 
the corresponding pdf attains a (q = 1.6)- Gaussian form for N < 2 16 (see Fig. HTk). No transition to a different 
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FIG. 8: (Color online) Plots of the ^-logarithm (inverse function of the ^-exponential ([7j)) vs. (z/a) 2 applied to our data of 
the normalized pdf of the (e = 0.9, \i — 1.6)-MacMillan map. N is the number of iterates, starting from Ni c randomly chosen 
initial condition in a square (0, 10 -6 ) x (0, 10 -6 ). For g-Gaussians this graph is a straight line, whose slope is —f3) for the right 
value of q. Note that the pdfs approach a true Gaussian (with (3 — 1) since q tends to 1 as N increases. 
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FIG. 9: Detailed evolution of the pdfs of the MacMillan map for e — 0.9, \i — 1.6, as N increases from 2 12 to 2 26 , respectively. 



QSS is detected, until the orbits diffuse to a wider chaotic region in phase space, for N < 2 18 . Let us observe in 
Fig. IU the corresponding pdfs of the rescaled sums of iterates, where even the tail of the pdf appears to converge to 
a (q = 1.6)-Gaussian (Fig. [TTb). For larger TV, further diffusion ceases as orbits "stick" to the outer islands, where 
the motion stays from there on. This only affects the tail of the distribution, which now further converges to a true 
(q = 1.6)- Gaussian representing this QSS up to N = 2 20 ). 

The remaining cases of Figures [3] and 2] can be viewed from a similar perspective. Indeed, the above analysis of the 
e = 1.2 example can serve as a guide for the (e = 0.5, /i = 1.6)- and (e = 1.1, /i = 1.6)-MacMillan maps, as well. In 
every case, the smallness of the L max but also the details of the diffusion process seem to play a key role in explaining 
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FIG. 10: Structure of phase space plots of the MacMillan map for parameter values e = 1.2 and \i — 1.6, starting from a 
randomly chosen initial condition in a square (0, 10 -6 ) x (0, 10 -6 ), and for N iterates. 




FIG. 11: (Color online) Pdfs of the rescaled sums of iterates of the MacMillan map for e = 1.2 and fi = 1.6 are seen to converge 
to a (q = 1.6)-Gaussian. This is shown in the panel (a) for the central part of the pdf (for N < 2 18 ) and in the panel (b) for the 
tail part (N > 2 18 ). Ni c is the number of initial conditions that have been randomly chosen from a square (0, 10~ 6 ) x (0,10~ 6 ) 



the convergence of pdfs to a g-Gaussian. What differs is the particular phase space picture that emerges and the 
number of iterations required to achieve the corresponding QSS. 

We conclude, therefore, that the dynamics of the MacMillan map for \i = 1.6 and e = 0.2,0.9, 1.8, where chaotic 
orbits evolve around the two islands of a single 'figure eight' chaotic region possess pdfs which pass rather quickly 
from a g-Gaussian shape to exponential to Gaussian. By contrast, the cases with e = 0.5, 1.1, 1.2 possess a chaotic 
domain that is considerably more convoluted around many more islands and hence apparently richer in "stickiness" 
phenomena. This higher complexity of the dynamics may very well be the reason why these latter examples have 
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motion initial conditions yo. In all 


cases, q x — 0.21, q y — 0.24, xo 


= -0.0049, xi = -0.5329, and yi = 0. 




yo 


ymax 


I 


0.00001 


0.00002 


II 


0.0001 


0.0003 


III 


0.001 


0.004 


IV 


0.01 


0.015 



QSS with g-Gaussian-like distributions that persist for very long. Even though we are not at an 'edge of chaos' 
regime where L max = 0, we suggest that it is the detailed structure of chaotic regions, with their network of islands 
and invariant sets of cantori, which is responsible for obtaining QSS with long-lived g-Gaussian distributions in these 
systems. 

III. FOUR-DIMENSIONAL CONSERVATIVE MAPS 

We now briefly discuss some preliminary results on the occurrence of QSS and nonextensive statistics in a 4- 
dimensional symplectic mapping model of accelerator dynamics [9]. This model describes the effects of sextupole 
nonlinearities on a hadron beam passing through a cell composed of a dipole and two quadrupole magnets that 
focuses the particles' motion in the horizontal (x)- and vertical (y)-directions |l0j. After some appropriate scaling, 
the equations of the mapping are written as follows: 

^n+l = 2c x X n X n —i P% n H~ y n (10) 

y n+ i = 2c y y n - y n -i + 2x n y n 

where p = P x s x /f3 y s y , c x , y = cos (2irq Xiy ) and s x , v = sin (2irq x ^y), q x ^ y is the so-called betatron frequencies and /3 Xi y are 
the betatron functions of the accelerator. As in [9], we assume that f3 XyV are constant and equal to their mean values, 
i.e. proportional to q~ y (q x = 0.21, q y = 0.24) and place our initial conditions at a particular point in 4-dimensional 
space associated with weak diffusion phenomena in the ^/-direction. In particular, our (xo,xi) — (—0.0049, —0.5329) 
coordinates are located within a thin chaotic layer surrounding the islands around a 5-order resonance in the x ni x n -\ 
plane of a purely horizontal beam, i.e. when y n = y n -i — 0. We then place our initial 2/1,2/0 coordinates very close to 
zero and observe the evolution of the y n s indicating the growth of the beam in the vertical direction as the number 
of iterations N grows. 

Let us observe this evolution in Fig. [12] separately in the x n +i,x n (first column) and y n +i,y n (second column) 
2-dimensional projections of our chaotic orbits. Clearly the behavior of these projections is very different: In the 
x-plane the motion keeps evolving in a thin chaotic layer around five islands, "feeding" as it were the (y n ,y n +i) 
oscillations, which show an evidently slow diffusive growth of their amplitude outward. 

In Table EH we list, for different initial values of 1/0 (yi — 0), the maximum amplitude of the y— oscilaltions, ymax-, 
while Fig. [T3l shows the corresponding pdfs of the normalized sums of iterates of the ?/n-variable. Note that, just as 
in the case of 2-dimensional maps, these distributions are initially q- Gaussian-\ike evolving slowly into triangular-\ike 
distributions, which finally approach Gaussians. In Fig. [13] we follow this evolution by performing four computations 
of TV = 2 19 iterates starting with yo which increases every time by a factor of 10. 

This similarity with the 2-dimensional case makes us suspect that the orbits of our 4-dimensional map also follow a 
sequence of weakly chaotic QSS, whose time-evolving features are depicted in plots of the ^-motion in Fig [T2l (second 
column), for increasing N. Note, for example, that one such QSS with a maximum amplitude of about 0.00001 is 
observed up to N > N = 2 19 , diffusing slowly in the ^/-direction. The pdf of this QSS is shown in the upper left 
panel of Fig [13] and has the shape of a g-Gaussian up to this value of N. However, for higher values of the yo initial 
condition, due to the sudden increase of the y n amplitudes at N = 2 20 , the "legs" of the pdf are lifted upward and 
the distribution assumes a more triangular shape. 

This rise of the pdf "legs" to a triangular shape is shown in more detail in Fig. [H] for initial conditions yo = 
10 -5 , 10 -4 , as the number of iterations grows to N = 2 20 . Clearly, the closer we start to yo = y\ = the more our 
pdf resembles a g-Gaussian, while as we move further out in the ^-direction our pdfs tend more quickly towards a 
Gaussian-like shape. This sequence of distributions is reminiscent of what we found for the 2-dimensional MacMillan 
map at (e = 0.9, /i = 1.6). Recall that, in that case also, a steady slow diffusion was observed radially outward, similar 



12 




-0.8 -0.4 X n 0.4 0.8 .3 1Q -4 y n ' ' 3.10" 4 



FIG. 12: The x n ,x n +i (first column of panels) and y n ,y n +i (second column of panels) projections of a chaotic orbit of (|10p , 
with q x = 0.21, q y = 0.24, x = -0.0049 and initial conditions xi = -0.5329, y = 0.0001 and yi = (case II of Table III). N 
represents the number of plotted iterates. 

to what was observed for the 4-dimensional map ([TO]) , which does not appear to be limited by a closed invariant curve 
in the x n , y n plane. 

One might wonder if it is possible to obtain for the 4-dimensional map (fTOj) also long-lived g-Gaussian pdfs of the 
type we found in the 2-dimensional MacMillan map. The likelihood of this occurrence is small, however, as all orbits 
we computed for the accelerator map (jTUJ) eventually escaped to infinity \ This implies that stickiness phenomena 
on island boundaries and sets of cantori are more dominant and tend to slow down diffusion more in the plane of 
2-dimensional maps like the MacMllan map than the 4-dimensional space of the accelerator map. It would, therefore, 
be very interesting to study, in a future paper, higher-dimensional maps whose chaotic orbits never escape to infinity 
(e.g. coupled standard maps) and compare their statistics with what we have discovered for the examples treated in 
the present paper. 

IV. CONCLUSIONS 

Our work serves to connect different types of statistical distributions of chaotic orbits (in the context of the Central 
Limit Theorem) with different aspects of dynamics in the phase space of conservative systems. What we have found, in 
several examples of the McMillan and Ikeda 2-dimensional area preserving maps as well as one case of a 4-dimensional 
symplectic accelerator map, is that g-Gaussians approximate well quasi-stationary states (QSS), which are surprisingly 
long-lived, especially when the orbits evolve in complicated chaotic domains surrounding many islands. This may be 
attributed to the fact that the maximal Lyapunov exponent in these regions is small and the dynamics occurs close 
to the so-called "edge of chaos" where stickiness effects are important near the boundaries of these islands. 

On the other hand, in simpler-looking chaotic domains (surrounding e.g. only two major islands) the observed QSS 
passes, as time evolves, from a g-Gaussian to an exponential pdf and may in fact become Gaussian, as the number 
of iterations becomes arbitrarily large. Even in these cases, however, the successive QSS are particularly long-lasting, 
so that the Gaussians expected from uniformly ergodic motion are practically unobservable. 

Interestingly enough, similar results have been obtained in N-dimensional Hamiltonian systems 0, [l3| describing 
FPU particle chains near nonlinear normal modes which have just turned unstable as the total energy is increased. 
Since these models evolve in a multi-dimensional phase space, the g-Gaussian pdfs last for times typically of the order 
10 6 , then pass quickly through the triangular stage and converge to Gaussians, as chaotic orbits move away from thin 



13 




! Case IV 




: y =o.oi 










z/c 



10 



FIG. 13: Pdfs of the normalized sums of iterates of the y-chaotic orbit of the 4-dimensional map (|10f) . for different yo. In all 
cases, q x — 0.21, q y = 0.24, xo = —0.0049, x\ — —0.5329 and y\ — 0. The number of (summed) iterates is N = 2 19 , and the 
number of randomly chosen initial conditions within an interval [O.9i/o,2/o] is Ni c = 10 5 . 




FIG. 14: (Color online) Pdfs of the normalized sums of iterates of the ^/-chaotic orbit of the 4-dimensional map, for different 
initial conditions yo and numbers of (summed) iterates N. Ni C is the number of randomly chosen initial conditions from an 
interval [0.9yo, yo]. In all cases, q x = 0.21, q y = 0.24, x = -0.0049, xi = -0.5329, and yi = 0. 
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layers to wider "seas" , where Lyapunov exponents are much larger. However, as long as the motion evolves near an 
"edge of chaos" region the distributions are g-shaped for long times, exactly as we found in the present paper. 

These conclusions are closely related to results obtained by other authors 0, who also study QSS occurring in 
low-dimensional Hamiltonian systems like 2-D and 4-D maps, but not from the viewpoint of sum distributions. They 
define a variance of momentum distributions representing a temperature-like quantity T(t) and show numerically that 
T(t) follows a "sigmoid" curve starting from small values and converging to a final value, which they identify as the 
Boltzmann Gibbs (BG) state. Although their initial conditions are spread over a wide domain and do not start from 
a precise location in phase space as in our studies, they also discover many examples of QSS which remain at the 
initial temperature for very long times, before finally converging to the BG state. 
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